#include <iostream>
#include <fstream>
#include <vector>
#include <unistd.h>
#include <iostream>

#include "MACRO.H"
#include "Particle.h"
#include "Vector.h"
#include "Color.h"
#include "ExpDisk.h"
#include "TreeForceFinder.h"


class Cell {
public:
	
	double m;
	Vector<double> R0;
	Vector<double> C;
	Vector<double> R1;
	Vector<double> A;
	
	vector<Particle*> p;
	int id;

	Cell():id(-2) { }
};


class FlatGrid {
public:
	
	Vector<int> n;
	Vector<double> R0, R1;
	
	TreeForceFinder tff;
	
	vector<Cell> cell;
	
	FlatGrid():
	R0( Vector<double>(-32.4,-18,0) ),
	R1( Vector<double>(+32.4,+18,0) ),
 	n( Vector<int>(360,200,1) )
	{ tff.eps = 0.08; }
	
	void create() {
		
		cell.resize((n.x+1)*(n.y+1)*n.z);
		
		Vector<double> DR((R1.x-R0.x)/n.x, (R1.y-R0.y)/n.y, 0);
		
		for (int j=0; j<=n.y; j++)
		for (int i=0; i<=n.x; i++) {
			
			Cell &cellijk = cell[i + (n.x+1)*j];
			cellijk.R0 = R0+Vector<double>(i*DR.x, j*DR.y, 0);
			cellijk.C  = R0+Vector<double>((i+0.5)*DR.x, (j+0.5)*DR.y, 0);
			cellijk.R1 = R0+Vector<double>((i+1)*DR.x, (j+1)*DR.y, 0);
			
		}
	}

	
	void get_acc(Particle &p) {
	
		Vector<double> DR(R1.x-R0.x, R1.y-R0.y, 0);
		Vector<double> iDR( DR.x!=0 ? n.x/DR.x : 0, DR.y!=0 ? n.y/DR.y : 0,  0 );
			
		Vector<double> R = p.R - R0;
			
		double dix = R.x*iDR.x;
		double diy = R.y*iDR.y;
		
		int ix = dix;
		int iy = diy;
		
		Vector<double> O1(dix-ix, diy-iy, 0);
		Vector<double> O0(1-O1.x, 1-O1.y, 0);
		
		
		Cell &cell00 = cell[ix+0 + (n.x+1)*(iy+0)];
		Cell &cell10 = cell[ix+1 + (n.x+1)*(iy+0)];
		Cell &cell01 = cell[ix+0 + (n.x+1)*(iy+1)];
		Cell &cell11 = cell[ix+1 + (n.x+1)*(iy+1)];
				
		double w[2][2];
				
		w[0][0] = O0.x*O0.y;
		w[0][1] = O0.x*O1.y;
		w[1][0] = O1.x*O0.y;
		w[1][1] = O1.x*O1.y;
				
		p.Ar.x =	cell00.A.x*w[0][0] +
					cell01.A.x*w[0][1] +
					cell10.A.x*w[1][0] +
					cell11.A.x*w[1][1];
					
		p.Ar.y =	cell00.A.y*w[0][0] +
					cell01.A.y*w[0][1] +
					cell10.A.y*w[1][0] +
					cell11.A.y*w[1][1];
											
		p.Ar.z	= 0;
	
	}
	
	
	void set_accs(vector<Particle> &p) {
		
		tff.cell0.R0 = 1.1*R0;
		tff.cell0.S = 1.1*R1-1.1*R0;
		
		tff.build_tree(p);
		//tff.cell0.print();
		
		//tff.test_tree(p);
		
		for (int i=0; i<cell.size(); i++)
			tff.calc_acc(cell[i].R0, cell[i].A);
		
	}
	

	void print() {
		printf("\n");
		printf(_B"FlatGrid::print()\n"B_);
		printf(" ├ cell.size() = %i\n", cell.size() );
		printf(" ├ cellsize = {%20.16e, %20.16e, %g}\n", (R1.x - R0.x)/n.x, (R1.y - R0.y)/n.y, 0);
		printf(" ├ R0 = {% 20.16e, % 20.16e, %g}\n", R0.x, R0.y, R0.z);
		printf(" └ R1 = {% 20.16e, % 20.16e, %g}\n", R1.x, R1.y, R1.z);
		printf("\n");
	}
	
	
};

class Regularizator {
	public:
	
	
	double stepsize;
	int nsteps;
	
	
	FlatGrid flatgrid;
	
	TreeForceFinder tff;
	
	Regularizator():
	stepsize(0.001),
	nsteps(2000),
	gp(popen("gnuplot --persist", "w"))
	{
		tff.eps = 0.01;
	}
		
	void regularize(vector<Particle> &p, vector<Particle> &p0) {
	
		printf("Create lookup acc grid ... "); fflush(stdout);
		
		flatgrid.create();
		flatgrid.set_accs(p0);
		
		printf("done.\n"); fflush(stdout);

		
		for (int frame=0, step=0; frame<340; frame++) {
		
			// full accs
			tff.build_tree(p);
			tff.calc_accs(p);
	
			// regular accs
			for (int i=0; i<p.size(); i++)
				flatgrid.get_acc( p[i] );
			
			
			if (30<frame && frame<310)
			for (int i=0; i<2; i++, step++ ) {

				// here we are suppressing irregular accs
				for (int i=0; i<p.size(); i++) {
					Vector<double> DA =	p[i].A - p[i].Ar; DA.z = 0;
					p[i].R-=stepsize*DA.unit();
				}
			
			}
			
			char fname[256];
			sprintf(fname, "/home/yurinds/tmp/voronoi/%04d.png", frame);
			plot(p, p0, step, fname);
			
			printf("%i %i\n", frame, step);
			
			
		}
		
	}
	
	
	void plot(vector<Particle> &p, vector<Particle> &p0, int step, const char *fname) {
		
		//double W = 1920;
		//double H = 1080;
		
		double W = 1700;//1920;
		double H = 956;//1080;
		
		double rm = 10;
		double lm = 10;
		double tm = 10;
		double bm = 10;
		
		double hr = 40;
		double wr = hr*(W-lm-rm)/(H-tm-bm);
		
		double boxsize = 0.15*(H-tm-bm)/flatgrid.n.y; // 1 - 0.165 2 - 0.15

		//gp = stdout;
		

		
		fprintf(gp, "set term wxt size %g, %g\n", W, H);
		
		//fprintf(gp, "set term pngcairo size %g, %g enhanced\n", W, H);
		//fprintf(gp, "set output '%s'\n", fname);
		
		
		fprintf(gp, "set xrange [0.5*%g:0.5*%g]\n", flatgrid.R0.x, flatgrid.R1.x);
		fprintf(gp, "set yrange [0.5*%g:0.5*%g]\n", flatgrid.R0.y, flatgrid.R1.y);
		
		fprintf(gp, "set lmargin at screen %g\n", lm/W) ;
		fprintf(gp, "set rmargin at screen %g\n", (W-rm)/W);
		
		//fprintf(gp, "set tmargin at screen %g\n", tm/H);
		//fprintf(gp, "set bmargin at screen %g\n", (H-tm)/H);
		
		//fprintf(gp, "set border lc -1 lw 2\n");
		
		fprintf(gp, "set xtics %g, %g\n", flatgrid.R0.x, (flatgrid.R1.x-flatgrid.R0.x));
		fprintf(gp, "set ytics %g, %g\n", flatgrid.R0.y, (flatgrid.R1.y-flatgrid.R0.y));
		
		fprintf(gp, "set mxtics %i\n", flatgrid.n.x);
		fprintf(gp, "set mytics %i\n", flatgrid.n.y);
		
		fprintf(gp, "set format x ''\n");
		fprintf(gp, "set format y ''\n");
		//fprintf(gp, "set object 1 rectangle from screen 0,0 to screen 1,1 fillcolor rgb '#FF5555' back\n");
		
		fprintf(gp, "set grid mxtic mytics back lt 1, lt 3 lc rgb '#EEEEEE'\n");
		
		//fprintf(gp, "set grid\n");
		
		fprintf(gp, "set size ratio -1\n");
		
		double s(0);
		for (int i=0; i<p.size(); i++)
			if (0<p[i].Ar.norm()) s += (p[i].A-p[i].Ar).norm()/p[i].Ar.norm();
		s/=p.size();	
		
		//fprintf(gp, "set label 1 'irregular accelerations' at screen 0.45, 0.35 tc lt 1 font 'Arial, 20' front\n");
		 
		fprintf(gp, "set label 100000 \"");
		fprintf(gp, "p.size() = %i\\n", p.size() ); //(number of points)
		fprintf(gp, "p0.size() = %i\\n", p0.size() ); //  (number of acc field generation points)
		fprintf(gp, "flatgrid.cell.size() = %i\\n", flatgrid.cell.size() ); //  (number of cells in acc lookup grid)
		fprintf(gp, "flatgrid.tff.eps = %g\\n", flatgrid.tff.eps );
		fprintf(gp, "regularizator.tff.eps = %g\\n", tff.eps );
		fprintf(gp, "regularizator.stepsize = %g\\n", stepsize );
		fprintf(gp, "s =%7.3f\\n", s );
		fprintf(gp, "step = %i\\n", step );
		fprintf(gp, "\" at screen 0.025, 0.93 front font 'Arial, 18'\n");
		

		fprintf(gp, "plot\ ");
		//fprintf(gp, "'-' u 1:2:3 w p pt 5 ps %g lc rgb var t '',\ ", boxsize);
		fprintf(gp, "'-' u 1:2:3 w p pt 7 ps 0.7 lc rgb var t '',\ ");
		//fprintf(gp, "'-' u 1:2:3:4:5 w vectors head filled lt 1 lw 2 lc rgb var t '',\ ");
		fprintf(gp, "'-' u 1:2:3:4 w vectors head filled lt 1 t''\n");

		/*
		// all grid cells
		for (int i=0; i<cell.size(); i++) {
			Color color; color.HSV2RGB(0.4, 0.10, 1.0);
			fprintf(gp, "%e %e %i\n", cell[i].C.x, cell[i].C.y, color.getColor() );
		}
		fprintf(gp, "e\n");
		*/
		
		// mesh generation points
		for (int i=0; i<p.size(); i++) {
			Color color; color.HSV2RGB(0.65, 1.0, 1.0);
			fprintf(gp, "%e %e %i\n", p[i].R.x, p[i].R.y, color.getColor());
		}
		fprintf(gp, "e\n");
		
		
		double f = 0.0001;
		
		/*
		// mesh generation vectors
		for (int i=0; i<p.size(); i++) {
			Color color; color.HSV2RGB(0.95, 1.0, 1.0);
			fprintf(gp, "%e %e %e %e %i\n", p[i].R.x, p[i].R.y, f*(p[i].A.x-p[i].Ar.x), f*(p[i].A.y-p[i].Ar.y), color.getColor());
		}
		fprintf(gp, "e\n");
		*/
		// acc vectors
		for (int i=0; i<1/*cell.size()/**/; i++) {
			Color color; color.HSV2RGB(0.0, 1.0, 1.0);
			fprintf(	gp, "%e %e %e %e\n", 
						flatgrid.cell[i].R0.x, flatgrid.cell[i].R0.y, 
						f*flatgrid.cell[i].A.x, f*flatgrid.cell[i].A.y );
		}
		fprintf(gp, "e\n");
		
		
	}
	
	private:	
	FILE *gp;
	
};



int main(int argc, char *argv[]) {
		
	srand48(50);
	
	printf("\n");
	
	printf("Create field particle distribution (p0): "); fflush(stdout);
	
	vector<Particle> p0(2000000);	
	ExpDisk expdisk;
	expdisk.rmax = 18;
	expdisk.set_rm(p0);
	
	for (int i=0; i<p0.size(); i++)  { 
		//p0[i].R.y = p0[i].R.z; 
		p0[i].R.z = 0; 
	}
	
	printf("Create particle distribution (p): "); fflush(stdout);
	
	vector<Particle> p(500);
	expdisk.set_rm(p);
	for (int i=0; i<p.size(); i++) { 
		//p[i].R.y = p[i].R.z;  
		p[i].R.z = 0; 
	}
	
	printf("Start regularization:\n"); fflush(stdout);
	
	Regularizator regularizator;
	regularizator.regularize(p, p0);
	
	return 0;
	
}